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Abstract. The effect of the intersite and interplane Coulomb interactions between the Dirac fermions on the 
formation of the Kohn-Luttinger superconductivity in bilayer doped graphene is studied disregarding the 
effects of the van der Waals potential of the substrate and both magnetic and non-magnetic impurities. The 
phase diagram determining the boundaries of superconductive domains with different types of symmetry of 
the order parameter is built using the extended Hubbard model in the Born weak-coupling approximation 
with allowance for the intratomic, interatomic, and interlayer Coulomb interactions between electrons. It 
is shown that the Kohn-Luttinger polarization contributions up to the second order of perturbation theory 
in the Coulomb interaction inclusively and an account for the long-range intraplane Coulomb interactions 
significantly affect the competition between the superconducting f— , p + ip— , and d + id—wa,ve pairings. It 
is demonstrated that the account for the interplane Coulomb interaction enhances the critical temperature 
of the transition to the superconducting phase. 

PACS. 74.20.Mn Superconductivity: nonconventional mechanisms - 74.25.Dw Superconductivity: phase 
diagrams - 74.78.Fk Superconducting multilayers - 81.05.ue Carbon-based materials: graphene 


1 Introduction 

In recent years, there has been an increased interest in the 
possibility of the development of the Cooper instability 
in graphene under appropriate experimental conditions. 
Although so far this possibility has not been confirmed, 
it was experimentally shown |TJ[2[3j[4l[5j[6] that graphene 
becomes superconducting when it is in a contact with or¬ 
dinary superconductors. This fact stimulated theoretical 
studies on possible implementation of the superconduct¬ 
ing phase in an idealized monolayer and bilayer graphene 
where the authors did not take into account the effect of 
nonmagnetic impurities and van der Waals potential of 
the substrate. 

Along with the numerous studies of this problem us¬ 
ing the electron-phonon mechanism [7l[8lf9l ll01lllj . pair¬ 
ing mechanisms caused by electron correlations mmm 
ns], and other exotic superconductivity mechanisms m 
Hz], some authors widely discuss the possibility of the de¬ 
velopment of Cooper instability in the above-mentioned 
systems using the Kohn-Luttinger mechanism which 
suggests the emergence of superconducting pairing in the 
systems with the purely repulsive interaction mmm 

mm- 
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As it was shown in [!M], the Cooper instability can 
occur in an idealized graphene single layer due to the 
strong anisotropy of the Fermi contour for Van Hove filling 
nvHi which, in fact, originates from the Kohn-Luttinger 
mechanism. According to the results obtained in m, this 
Cooper instability in graphene evolves predominantly in 
the d— wave channel and can be responsible for the critical 
superconducting transition temperatures up to T c ~ 10 1\ , 
depending on the proximity of the chemical potential level 
to the Van Hove singularity. The theoretical analysis of 
the competition between the ferromagnetic and supercon¬ 
ducting instabilities showed [25| that the tendency to su¬ 
perconductivity due to strong modulation of the effective 
interaction along the Fermi contour, i.e., due to electron- 
electron interactions alone, prevails. In this case, the su¬ 
perconducting instability evolves predominantly in the /- 
wave channel. 

The competition between the Kohn-Luttinger super¬ 
conducting phase and the spin density wave phase at the 
Van Hove Filing and near it in the graphene single layer 
was analyzed in mm using the functional renormaliza¬ 
tion group method. It was found that superconductivity 
with the d + id— wave symmetry of the order parameter 
prevails in a large domain near the Van Hove singularity, 
and a change in the calculated parameters may lead to 
a transition to the phase of the spin density wave. Ac- 
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cording to E2, far away from the Van Hove singularity, 
the long-range Coulomb interactions change the form of 
the d + id —wave function of a Cooper pair and can facili¬ 
tate superconductivity with the /—wave symmetry of the 
order parameter. The competition between the supercon¬ 
ducting phases with different symmetry types in the wide 
electron density range 1 < n < nvH in the graphene sin¬ 
gle layer was studied in [281129] . It was demonstrated that 
at intermediate electron densities the long-range Coulomb 
interactions facilitate implementation of superconductiv¬ 
ity with the /—wave symmetry of the order parameter, 
while at approaching the Van Hove singularity, the su¬ 
perconducting pairing with the d + id—symmetry type 
evolves mm- 

The conditions for the Kohn-Luttinger superconduct¬ 
ing pairing was analyzed also in graphene bilayer 15011511 
mm- According to the results of [3J], the ferromag¬ 
netic instability near the Van Hove singularities domi¬ 
nates over the Kohn-Luttinger pairing in graphene bilayer. 
It should be noted, however, that in these calculations 
only the Coulomb repulsion of electrons on one site was 
taken into account. Authors of [35] calculated the screen¬ 
ing function of Coulomb interaction in the doped and un¬ 
doped bilayer graphene in the random phase approxima¬ 
tion (RPA). They established that the static polarization 
operator in the doped regime contains the singular part 
(the Kohn anomaly) that significantly exceeds one calcu¬ 
lated for monolayer or 2D electron gas. As it is known, 
the Kohn anomaly EM57J facilitates the effective attrac¬ 
tion between two particles, inducing a contribution that 
always exceeds the repulsive contribution connected with 
the regular part of the polarization operator for the angu¬ 
lar momenta l ^ 0 oi two particles [18]. Therefore, one can 
expect that the critical superconducting temperature T c 
in an idealized bilayer can exceed the corresponding value 
for graphene monolayer. 

Additionally, it was shown in papers [38li39] that the 
value of T c can be increased in the framework of the Kohn- 
Luttinger mechanism even for low carrier densities if the 
spin-polarized two-band situation or a multilayer system 
is considered. In this situation, the role of the pairing spins 
”up” is played by electrons of one band (layer), while the 
role of the screening spins ” down” is played by electrons of 
another band (layer). Coupling between the electrons from 
the two bands occurs owing to the interband (interlayer) 
Coulomb interaction. In this case, the following mecha¬ 
nism is possible: electrons of one sort form a Cooper pair 
by polarizing the electrons of another sort [381139] . This 
mechanism can be realized also in quasi-2D systems. 

In this paper, in the Born weak-coupling approxima¬ 
tion, we consider the Kohn-Luttinger superconducting pair¬ 
ing in an idealized graphene bilayer. We calculate the 
phase diagram, which reflects the competition between 
the superconducting phases with different types of the 
symmetry of the order parameter, taking into account the 
second-order contributions in the Coulomb interaction to 
the effective interaction of electrons in the Cooper chan¬ 
nel. We analyze modification of the phase diagram with 
allowance for the Coulomb repulsion between electrons of 



Fig. 1 . Crystalline structure of bilayer graphene. Atoms Ai 
and Bi in the lower monolayer are shown by red and black 
circles; black and green circles in the upper layer correspond 
to atoms A 2 and B 2 . Intraplane electron hoppings are labeled 
by fi and £ 2 ; 71 , 73 and 74 show the interlayer hoppings. 


the same, of the nearest, and of the next-to-nearest carbon 
atoms in a single layer, as well as the interlayer Coulomb 
interactions. We demonstrate the importance of taking 
into account the Coulomb repulsion of electrons on dif¬ 
ferent crystal lattice sites and in different layers of bilayer 
graphene. The account of Coulomb repulsion changes the 
phase diagram of the superconducting state and, under 
certain conditions, increases the critical temperature. 


2 Theoretical model 

We consider an idealized graphene bilayer, assuming that 
two layers are arranged in accordance with the AB type, 
i.e., one layer is rotated on 60° relative to the other one [40] 
E]. Let us choose the arrangement of the sublattices in 
the layers in such a way that the sites from different layers 
located one above another belong to the sublattices Ai 
and A -2 respectively, while the other sites belong to the 
sublattices B\ and B 2 (Fig. [I])- In the Shubin-Vonsovsky 
(extended Hubbard) model [H], the Hamiltonian for the 
graphene bilayer which takes into account electron hop¬ 
pings between the nearest and next-to-nearest atoms, as 
well as the Coulomb repulsion between electrons of the 
same and of the adjacent atoms and the interlayer Coulomb 
interaction of electrons, in the Wannier representation has 
the form: 



tl 'y ^ ( a j f a biJ+S,(7 + (3'2/ct^ 2,/-(5,(T + h.c.) 
fSa 
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As a result, Hq acquires the form 

4 

H 0 = EE Eik(%\ kcr Oiikcr ■ (5) 

2=1 ka 

Since the results of ab initio calculations for graphite [33) 
03 showed a very small value of the interlayer hopping 
parameter 74 , hereinafter we assume that 74 = 0. Then, 
the four-band energy spectrum of the graphene bilayer is 
described by the expressions 


E lk = e ± \j'A k ± \Z~Bk — t 2 f kl ( 6 ) 

A k = — ^2 a 2 + 4|6fc | 2 + 2 |dfc| 2 ^, 

B k = ^{\d k \ 2 [\d k \ 2 ^ 2a 2 +4|6 fc | 2 ) + a 4 +4a 2 |6 fc | 2 

+4a6|d fc + 4afef4), 

a = 71 , b k = hu k , d k = 7 3 u k , 


where the following notation has been introduced: 

'3 


(y/l 

f k = 2 cos(V 3 k y ) + 4cos( 


cos 


( 7 ) 


In (0]) |3]), the operators a{j 0 .(ai/ 0 -) create (annihilate) an 
electron with the spin projection a = ± 1/2 at site / of the 
sublattice A±; h^ a = a\ fiy aif a denotes the operators of 
the numbers of fermions at the / site of the sublattice 
A\ (analogous notations are used for the sublattices A 2 , 
B 1 , and B 2 ). Vector 5(—5) connects the nearest atoms 
of the hexagonal lattice of the lower (upper) layer. In¬ 
dex i = 1,2 in Hamiltonian 0 denotes the number of 
layer. We assume that the one-site energies are identical 
{eAi = £Bi = e) and the position of the chemical poten¬ 
tial p and number of carriers n in graphene bilayer can 
be controlled by a gate electric field. In the Hamiltonian, 
t\ is the hopping integral between the neighboring atoms 
(hoppings between different sublattices), t 2 is the hopping 
integral between the next-to-nearest neighboring atoms 
(hoppings in the same sublattice), U is the parameter of 
Coulomb repulsion between electrons of the same atom 
with the opposite spin projections (Hubbard repulsion), 
and Vj and V 2 are the Coulomb interactions between elec¬ 
trons of the nearest and the next-to-nearest carbon atoms 
in a single layer. The symbol (( )) indicates that sum¬ 
mation is made only over next-to-nearest neighbors; the 
symbols 71 , 73 , 74 denote the parameters of the interlayer 
electron hoppings (Fig. [TJ , and Gi, G 3 and G 4 are the 
interlayer Coulomb interactions between electrons. 


Uk = J2 ^ = e ~ ikx + 2e ^ cos ( Y^) > ( 8 ) 

\u k \ = y/3 + f k . (9) 


In this paper, the conditions for the implementation 
of the Kohn-Luttinger superconductivity are analyzed by 
considering the situation when upon doping of the graphene 
bilayer the chemical potential falls into the two upper en¬ 
ergy bands E\ k and E 2k (Fig.0^). Then, if 74 7 ^ 0 and the 
inequality p > 71 is valid, the Fermi contour will consist 
of two lines (Fig. 0 }d) in the vicinity of each Dirac point 
for the electron densities 1 < n < nyH, where n is the 
electron density calculated per atoms of one layer. 

The coefficients of the Bogolyubov transformation can 
be found from the system of homogeneous equations 


/ Xi a b* k 0 \ 

a Xi 0 b k j 

b k 0 Xi d k I 

V 0 b k d k Xi ) 


/ w a \ 
Wi2 I 
Wi3 I 

\WilJ 


= 0 , 


( 10 ) 


where Xi = E ik - £ + t 2 f k . 

In the Bogolyubov representation, the Hamiltonian Hi nt 
© in terms of the operators ai ka , a 2ka , a 3kcr and a± k(J 
reads as follows: 

Hint = El Elj.im(k, p\q, s)al ka Oj p(J ai qa a msa 


We diagonalize the Hamiltonian Hq using the Bogolyubov 
transformation 


= wa(k)a lk v + w i2 (k)a 2ka 


ijlmcr 

hpqs 

x S(k + p — q — s) ( 11 ) 

+ Y E ^Im(k.,p\q,s)al kt al pi ai qi a mst 

ijlm 


+ w i3 (k)b lka + w i4 ,(k)b 2ka , i = 1,2, 3,4. 


x S(k +p — q- s), 


(Xikcr 


( 4 ) 
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a) 




Fig. 2. (a) Energy bands of graphene bilayer in the vicinity 
of the Dirac points and (b) formation of the Fermi contour at 
= 0, 71 = 0.12, 73 = 0.1, and [i = 0.7 (all the parameters 
are given in units of |ti|). 

where 5{x) is the Dirac delta-function and p\q, s) 

and r^. lm (k,p\q, s) are the initial amplitudes. The quan¬ 
tity 

( k ’P\<h S ) 

= ^(Vij;lm(k,p\q,s) + Vji- t ml{p, k \s,q ) 

+ G S;L( fc ’Pl9> s ) + G Smz(P’ k \ s . <?) 

+ G i%m( k ’P\<l’ S ) + G fi-,ml(Pi k \ S ’ ( l) 

+ G Sm( fe -Pl9- S ) + G Smz(P> «)) , ( 12 ) 

%Zm (fe,p|g,s) = 

+W^-pW.2 (fc)w^-4 (P)UZ* 4 (?) W* m2 (s)) 


+ -y X] f<t-P W ir{ k )Wjr{P) W lr{q)w* mr {s ), (13) 

r—1 

G gJ m (fc,p|g,s) = G iw a (feH- 2 (pM* 2 (g)w;^i(s), (14) 

G S;L ( fc ’Plg, S ) 

= G 3 Uq-pW i 3(fc)w j 4(p)w; 4 (g)u4 3 (s), (15) 

G ij;Zm ( fc >Pk> s ) = G 4(<_pWil(fe)w i 4(p)w; 4 (g)u4 1 (s) 

+u q - p w a {k)w j3 (p)w* a {q)w* rn2 {s) s j (16) 

corresponds to the intensity of the interaction of fermions 
with parallel spin projections, while the quantity 

r ij-im (k,p\q,s) = U iji i m (k,p\q,s) 

H“ Vij\lm (k,p\q,s) + Vji- m i{p,k\s,q) 

+ G S;L( fe ’Pl9’ S ) + G Smz(P> k \ S > 9) 

+ G <Sm P l 9 ’ S) + G S?mZ (P. <?) 

+ G iSrri( fc >Pl9> s ) + G 5ijLl(P. 9), (17) 

Uij-im {k,p\q,s) 

4 

= UJ2 Wir (k)wjr (p)w* r (q)w^ nr (s) (18) 

r=l 

describes the interaction of fermions with antiparallel spin 
projections. Indices correspond to the number of 

the energy band and acquire the values 1, 2, 3, or 4. 

3 Effective interaction and equation for the 
superconducting order parameter 

In this paper, we use the Born weak-coupling approxima¬ 
tion, in which the hierarchy of model parameters has the 
form 

W > U > y > V 2 > Gi > G 3 , G 4 , (19) 

where W is the bandwidth in graphene bilayer (0. In 
the calculation of the scattering amplitude in the Cooper 
channel, the condition dD allows us to limit the con¬ 
sideration to only the second-order diagrams in the ef¬ 
fective interaction of two electrons with opposite values 
of the momentum and spin and use the quantity T(p, k) 
for it. Figure [3] depicts the sum of diagrams which de¬ 
termines r(p,k). Here, solid lines correspond to Green’s 
functions for the electrons with opposite spin projections 
(light arrows) and — \ (black arrows). The first dia¬ 
gram describes the initial interaction of two electrons in 
the Cooper channel. Here, the wavy lines correspond to 
the initial interaction. The next four diagrams in Fig. [3] 
correspond to the second-order scattering processes Sr(p, k) 
and describe the polarization effects of the filled Fermi 
sphere. In the diagrams, the presence of solid lines with¬ 
out arrows means the summation over the spin projections 
values. 
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Fig. 3. The sum of diagrams for the effective interaction of 
electrons in the Cooper channel for graphene bilayer. Solid 
lines correspond to Green’s functions for electrons with spin 
projections + 1 and — | and energies corresponding to the en¬ 
ergy bands Ei, Ej, Ei, and A m . Subscripts i and j acquire the 
values of 1 or 2, whereas l and m acquire the values of 1, 2, 3, 
or 4. Momenta q\ and q 2 are defined by (EJ). 


The possibility of the Cooper pairing is determined by 
the features of the energy structure and the effective inter¬ 
action of electrons near the Fermi level [33] . If we assume 
that the chemical potential in doped graphene bilayer is 
located in the two upper bands Eik and E 2 k (Fig- E), 
we can consider the situation in which the initial and final 
momenta of electrons in the Cooper channel also belong 
to the two upper bands and analyze the conditions for the 
Kohn-Luttinger superconducting pairing. At that, indices 
i and j in the diagrams (Fig. [3]) will acquire the values of 
1 or 2 . 

Introducing the analytical expressions for the diagrams, 
we get the effective interaction in the form 


r(p,k) = r 0 {p,k) + sr(p,k), ( 20 ) 

r 0 (p,k) = r£. j:j (p,-p\-k,k), ( 21 ) 

5 f(p,k) = jj Y r U;jm(P,<li\ - fc,Pi) (22) 

Z,ra,pi 

2 _^ r 

l,m,pi 

X [ r il,mM^-P\P^~ k ) - r ii;jm(Q l.“Pl “ 


+ r iijrn(Ql,-P\ - k >Pl) 


r fm-AP^PA k ^qi) - r iL,lj(P’Pl\qi’ k ) } Xl,m(qii Pi)- 


Here, we use the notations for the generalized susceptibil¬ 
ities 


Xl,m{k,p) = 


f(Eik ) ~ f(E mp ) 
E m p Elk 


(23) 


where /( x) = (exp(^#) + 1) 1 is the Fermi-Dirac func¬ 
tion and the energies Eik are defined by the expressions ©■ 


Additionally, we have introduced the following notations 
for the combinations of the momenta 

qi=Pi+p-k, q 2 =pi-p k. (24) 

The renormalized expression for the effective interac¬ 
tion allows us to analyze the conditions for the occurrence 
of superconductivity in the system. It is known 1-1 bl that 
the development of the Cooper instability can be estab¬ 
lished from the consideration of the homogeneous part of 
the Bethe-Salpeter equation. At that, the dependence of 
the scattering amplitude r(p,k) on momentum k is fac¬ 
torized and we get the integral equation for the supercon¬ 
ducting order parameter A(p). After the integration over 
the isoenergetic contours, the problem of the Cooper in¬ 
stability can be reduced to the eigenvalue problem 

I57 | l3g p g ] 


I? / = (25) 

e q =p 

where the eigenvector is the superconducting order pa¬ 
rameter A(q) and the eigenvalues A satisfy the relation 
A -1 ~ ln(T c /IF). Here, the momenta p and q belong to 
the Fermi surface and vp{q ) is the Fermi velocity. Equa¬ 
tion El is solved in accordance with the common scheme 
described in The integration is fulfilled with the 

allowance for the fact that the Fermi contour near each 
Dirac point consists of two lines (Fig. [2}j). 


4 Results and discussion 

Let us consider the phase diagram of the superconduct¬ 
ing state of the graphene bilayer and the modifications of 
this diagram in the different regimes obtained by solving 
Eq. ®. When building the phase diagram, we divided 
the multisheet Fermi contour into 180 intervals and the 
Brillouin zone of the graphene bilayer, into 5 • 10 4 cells. 
It was established that the chosen method of division is 
sufficient for the correct description of the dependence of 
the effective coupling constant A on the electron density 
n J2F. Based on the obtained dependences A(n) for differ¬ 
ent values of the intersite V\ and interplane G\ , G 3 and G '4 
Coulomb interactions, we built the phase diagrams of the 
Shubin-Vonsovsky model for bilayer graphene, which re¬ 
flect the competition between the superconducting phases 
with different types of symmetry of the order parameter. 

So far, there has been no agreement regarding the val¬ 
ues of parameters of the intra- and interplanar Coulomb 
interactions in the graphene bilayer. The ab initio calcu¬ 
lations for graphite J5D] showed that the value of Hubbard 
repulsion is U = 8.0 eV, which is consistent with the es¬ 
timation made in m and contradicts the intuitively ex¬ 
pected small value of U and weak-coupling limit U < W 
(it is known [52] that t\ « 2.8eV). The authors of [50] 
calculated the parameters of Coulomb repulsion between 
electrons of the nearest and the next-to-nearest carbon 
atoms: V\ = 3.9 eV and V 2 = 2.4 eV, respectively. At the 
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t =1 t =0 v =0 y =0 U=2 V=0.6V G=G=G=0 

1 2 1 1 "3 2 113 4 



Fig. 4. Phase diagram of the superconducting state of the 
graphene bilayer shown as a function of the variables ”n — Vi” 
at t 2 = 0, 71 = 73 = 0, U = 2, Gi = G 3 = G 4 = 0 for (a) 
V 2 = 0 and (b) V 2 = O. 6 V 1 (all the parameters are in units of 
|£i |). For all the points on the same thin blue line, the value of 
|A| is constant and marked with the corresponding number. 


same time, the other authors (see, for example, [53]) con¬ 
sider these parameters to be much smaller. The authors 
of [15] mentioned that the estimation of the parameters of 
Coulomb interaction, including the Hubbard repulsion, in 
the graphene bilayer strongly depends on the calculation 
scheme which is used. In our calculation, we apply the pa¬ 
rameter hierarchy m, which allows us to use the Born 
weak-coupling approximation. For interlayer hopping pa¬ 
rameters 71 and 73 , we use the values similar to those 
determined in [43U44 for graphite. 

First, let us consider the limiting case when the bilayer 
energy spectrum is described by the only one hopping pa¬ 
rameter (ti 7 ^ 0, t 2 =71 =73 =0). The Hubbard repul¬ 
sion is also taken into account U = 2 (hereinafter, all the 
parameters are given in units of |ti |). The Coulomb repul¬ 
sion between electrons (Vi 7 ^ 0) of the neighboring carbon 
atoms in the same layer is taken into account as well. At 



Fig. 5. Modification of the dependence of the density of elec¬ 
tron states of the graphene bilayer from electron density cal¬ 
culated for atoms in the same layer with respect to electron 
hoppings to the next-to-nearest neighboring atoms. 


the same time, the interlayer Coulomb interactions are not 
taken into account (Gi = G 3 = G 4 = 0). Thus, in the cho¬ 
sen regime, the graphene bilayer consists of two isolated 
single layers. The phase diagram of the superconducting 
state shown as a function of the variables ”n — Vi” for 
this case is presented in Fig. HJi. It can be seen that the 
phase diagram comprises three regions. At low electron 
densities n, the ground state of the system corresponds to 
the superconductivity with the d + id— wave symmetry of 
the order parameter, which is described by the 2D repre¬ 
sentation E 2 , the contribution to which is determined by 
the harmonics 

9m +td H^ ) = —i= (A sin (2m + 2)<j) + B cos (2m + 2 )cj>), 

V k 

where subscripts m run over the values for which the co¬ 
efficients (2m + 2 ) are not multiples of 3. At the inter¬ 
mediate electron densities, the superconducting /—wave 
pairing is implemented, the contribution to which is de- 

/ / \ 1 

termined by the harmonics g™ vt>) = sin (6m + 3)</> 
(here m £ [ 0, 00 )), while the contribution of the harmonics 

( f \ 1 

9m (<P) = — 7 = cos ( 6 m + 3)/) is absent. At the large val- 
v 7r 

ues of n, the domain of the superconducting d + id— wave 
pairing occurs m • With the increase of the parameter V\ 
of the intersite Coulomb interaction, in the region of small 
values of n, the d+id— wave pairing is suppressed and the 
pairing with the /—wave symmetry of the order parame¬ 
ter is implemented. Thin blue lines in Fig. [I] are the lines 
of the equal values of the effective coupling constant |A|. 
It can be seen that in this case in the proximity of the 
Van Hove filling tivh (solid curve in Fig. [5]) the effective 
coupling constant attains the values |A| = 0.1. 

It should be noted that to avoid the summation of the 
parquet diagrams [M1I5M5H] . we do not analyze here the 
electron density regions that are very close to the Van 
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t=l t 2 =0 Yj=0.12 y 3 =0-1 
a) U =2 V n =0 G=0.5V 1 G 3 =G 4 =0.4V 1 



b) t=l t=0 Yj=0.12 y 3 =0.1 

U=2 Y=0.6V G 1 =0.5V [ G 3 =G 4 =0.4V 1 



Fig. 6. Phase diagram of the superconducting state of the 
graphene bilayer shown as a function of the variables ”n — Vi” 
for t 2 = 0, 7 i = 0 . 12 , 73 = 0.1, - 2 , G'i = 0.514, G 3 = 
G± = 0.4Vi, at (a) V 2 = 0 and (b) V 2 = 0.6Vi (all the param¬ 
eters are given in units of |fi|). Thin blue curves are the lines 
of the constant values of |A|. 


Hove singularity in the density of electron states of bi¬ 
layer graphene (Fig. [5j. For this reason, the boundaries 
between different domains of the implementation of the 
Kohn-Luttinger superconducting pairing, as well as the 
lines of the equal value of |A| that are very close to the 
Van Hove singularity are indicated in the phase diagram 
by the dashed lines. 

Thus, in the numerical calculation for the graphene 
bilayer for the chosen parameters, we made the limiting 
transition to the results obtained by us previously for the 
graphene monolayer mm- 

Let us consider the modification of the phase diagram 
for the isolated graphene single layers with regard to the 
long-range intraplane Coulomb interactions between elec¬ 
trons V 2 . It can be seen in Fig. [4 )d for the fixed ratio be¬ 
tween the parameters of the long-range Coulomb interac¬ 
tions V 2 = O. 6 V 1 that when V 2 is taken into account, the 


phase diagram changes qualitatively. This change involves 
the suppression of a large domain of the superconduct¬ 
ing state with the /—wave symmetry at the intermediate 
electron densities and the implementation of the super¬ 
conducting pairing with the p + ip— wave symmetry of 
the order parameter. In addition, when V 2 is taken into 
account, the effective coupling constant increases to the 
value |A| = 0.3. 

Now, let us consider the modification of the phase di¬ 
agram of the superconducting state with respect to the 
interplanar interactions. When the interlayer electron hop¬ 
pings 71 = 0.12 and 73 = 0.1 are taken into account while 
the other parameters being the same as in Fig. [2 the 
phase diagram of the graphene bilayer remains nearly un¬ 
changed. 

Inclusion of the Coulomb interaction G\ in the consid¬ 
eration weakly shifts the boundaries of the / 1 —wave and 
d + id— wave pairing in the phase diagram in Fig. [I] and 
does not affect the absolute values of A. Figure [G] shows the 
effect of taking into account the interlayer Coulomb inter¬ 
actions G 3 and G4. Figure IHJi shows the phase diagram of 
the Shubin-Vonsovsky model for the graphene bilayer for 
the set of parameters t - 2 = 0, 71 = 0.12, 73 = 0.1, U = 2 
and V 2 = 0 for the chosen ratios between the interlayer and 
intersite Coulomb interactions Gi = 0.5Vi, G3 = G4 = 
0.4Vi, according to the hierarchy of the parameters m- 
The calculation shows that the separate increase of the 
parameters G 3 and G4 suppresses the d + id—wave pair¬ 
ing and, at the same time, broadens the /—wave pairing 
region at small electron densities. The superconducting 
d+id— phase is suppressed the most effectively by enhanc¬ 
ing the parameter G4 of the interlayer Coulomb interac¬ 
tion. When the interactions G3 and G4 are simultaneously 
taken into account (Fig. [2^), then along with the intensive 
suppression of the superconducting d+id— wave pairing at 
small electron densities and the implementation of the su¬ 
perconductivity with the /—wave symmetry of the order 
parameter, the growth of the absolute values of effective 
coupling constant A is also observed. 

Figure [5Jd depicts the phase diagram of the graphene 
bilayer calculated for the same parameters as in Fig. [ 6 jr 
but with respect to the long-range intraplane Coulomb 
repulsion between electrons V 2 . Comparison of Figs. (5J) 
and Eh shows that the account for G3 ^ 0 and G4 / 0 
leads to the strong competition between the d + id— wave 
and p + ip—wave pairings with the significant suppression 
of the p+ip— wave pairing in the region of the intermediate 
electron densities. In this case, in the remained region of 
the p + ip— wave pairing, |A p +i P | slightly exceeds |A/|. 

The account for electron hoppings to the next-to-nearest 
carbon atoms t 2 does not qualitatively affects the compe¬ 
tition between the superconducting phases (Fig. [ 6 ]). Fig¬ 
ure Q] depicts the phase diagram of the graphene bilayer 
obtained for the parameters t 2 = 0 . 1 , 71 = 0 . 12 , 73 = 
0.1, U = 2, Gi = 0.514 and G 3 = G 4 = 0.414. Such a be¬ 
havior of the system is explained by the fact that switching 
on of the hoppings t - 2 > 0 or t 2 < 0 for the graphene bi¬ 
layer, similarly to the case of the monolayer investigated 
by us in [M 23 j, does not significantly modify the den- 
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Fig. 7. Phase diagram of the superconducting state of the 
graphene bilayer shown as a function of the variables ”n — Vi” 
at f 2 = 0.1, 7 i = 0 . 12 , 73 = 0.1, U = 2 , Gi = 0.5Vi, G 3 = 
G 4 = 0.4Vi (all the parameters are in units of |ti|). Thin blue 
curves are the lines of the constant values of |A|. 


sity of electron states in the carrier concentration regions 
between the Dirac point and both points uvh (Fig. [5]). 
However, it can be seen in Fig. [7] that the account for 
the hoppings f 2 leads to an increase of the effective in¬ 
teraction in the absolute values and, consequently, to the 
higher superconducting transition temperatures in an ide¬ 
alized graphene bilayer. 

It should be noted that the Kohn-Luttinger supercon¬ 
ductivity in the graphene single layer and bilayer never 
develops near the Dirac points. The calculations show that 
in the vicinity of these points, where the linear approxima¬ 
tion for the energy spectrum of the graphene single layer 
and the parabolic approximation for the spectrum of the 
graphene bilayer work pretty well, the density of states is 
very low and the effective coupling constant |A| < 10~ 2 . 
The higher values of |A|, which are indicative of the de¬ 
velopment of the Cooper instability, arise at the electron 
densities n > 1.15. However, at such densities, the energy 
spectrum of the bilayer along the direction KM of the 
Brillouin zone (Fig. [2jo) already significantly differs from 
the Dirac approximation. 

5 Conclusions 

In the work, we have analyzed the conditions for the Kohn- 
Luttinger superconductivity in a semimetal with the Dirac 
spectrum using as an example an idealized graphene bi¬ 
layer, disregarding the van der Waals potential of the sub¬ 
strate and both magnetic and non-magnetic impurities. 
The electronic structure of graphene bilayer is described in 
the Shubin-Vonsovsky model taking into account not only 
the Coulomb repulsion of electrons of the same carbon 
atom, but also the intersite and interlayer Coulomb inter¬ 
actions. It was shown that in such a system, the Kohn- 
Luttinger polarization contributions lead to the effective 


attraction between electrons in the Cooper channel. The 
constructed superconducting phase diagram of the system 
determines the Cooper pairing domains with the different 
types of the symmetry of the order parameter, depend¬ 
ing on the intersite Coulomb interactions and the elec¬ 
tron densities. The analysis of the phase diagram showed 
that the inclusion of the Kohn-Luttinger renormalizations 
up to the second order of perturbation theory inclusively 
and the allowance for the long-range Coulomb interac¬ 
tions Vi and V 2 determine, to a considerable extent, the 
competition between the superconducting phases with the 
/—wave, p + ip—wave, and d + id-wave types of the sym¬ 
metry of the order parameter. They also lead to a signifi¬ 
cant increase in the absolute values of the effective inter¬ 
action. It was shown that the allowance for the interlayer 
Coulomb interactions G 3 and G 4 , as well as for the distant 
electron hoppings f 2 , leads to an additional increase in the 
effective interaction and, hence, to the higher supercon¬ 
ducting transition temperatures in an idealized graphene 
bilayer. 

Our calculation showed that the Kohn-Luttinger mech¬ 
anism can lead to the superconducting transition temper¬ 
atures T c ~ 20 -V 40 K in an idealized graphene bilayer. 
Contrary to these rather optimistic estimations, in real 
graphene, as it was mentioned in Introduction, supercon¬ 
ductivity has not been found yet. This material is only 
close to superconductivity. 

For a few reasons, the results of the theoretical cal¬ 
culations reported here can differ from the experimental 
situation. First, we did not take into account the effect of 
the van der Waals potential of the substrate |!571f551f5?| . It 
seems that the effect of this potential should be weakened 
with the increase of number of layers. However, even in 
the multilayer systems the van der Waals forces can de¬ 
grade the conditions for the development of the Cooper 
instability. 

Second, as we mentioned in Section 4, there has been 
no agreement regarding the values of the parameters of 
the intraplane and interplanar Coulomb interactions in 
the graphene bilayer in the literature. In this work, we 
used the values of the intraplane Coulomb interactions 
that are close to those obtained from the ab initio calcu¬ 
lation in [50] for graphite. The values of the interplanar 
Coulomb interactions were chosen to satisfy the hierarchy 
of the parameters of the Born weak-coupling approxima¬ 
tion. 

Third, in our calculations, we considered a pure graphene 
bilayer with the ideal structure, whereas the real material 
contains numerous impurities and structural defects. It 
is well known that, in contrast to the traditional s —wave 
pairing, for the anomalous pairing with the /-wave, p+ip- 
wave, and d + id -wave symmetries of the order parameter, 
nonmagnetic impurities and structural defects can destroy 
the superconducting order [60]. 

In addition, we should mention one more possible rea¬ 
son for the discrepancy between the theoretical calcula¬ 
tions on superconductivity in graphene and the experi¬ 
mentally observed situation. In recent paper m, the ef¬ 
fect of quantum fluctuations (T = 0) on the graphene 
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layers was investigated. It was shown that these fluctua¬ 
tions initiate the logarithmic corrections to the moduli of 
elasticity and bending of the layers. In other words, ac¬ 
cording to EH, the quantum fluctuations connected with 
the bending vibrations of the graphene layers can lead to 
the situation when the electrons do not move along the 
atomically smooth layers but along the strongly curved 
string-like trajectories, as in quantum chromodynamics. 
This situation requires further investigations, although in 
this case the superconductivity is not at all excluded and 
even can be enhanced by the exchange of bending vibra¬ 
tion quanta between the pairing electrons. 
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